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Abstract. Mechanical and sliding instabilities are the two processes which 
may lead to breaking off events of large ice masses. Mechanical instabilities 
mainly affect unbalanced cold hanging glaciers. For the latter case, a predic- 
tion could be achieved based on data of surface velocities and seismic activity. 
The case of sliding instabilities is more problematic. This phenomenon occurs 
on temperate glacier tongues. Such instabilities are strongly affected by the 
subglacial hydrology: melt water may cause (i) a lubrication of the bed and 
(ii) a decrease of the effective pressure and consequently a decrease of basal 
friction. Available data from Allalingletscher (Valais) indicate that the glacier 
tongue experienced an active phase during 2-3 weeks with enhanced basal mo- 
tion in late summer in most years. In order to scrutinize in more detail the 
processes governing the sliding instabilities, a numerical model developed to 
investigate gravitational instabilities in heterogeneous media was applied to 
Allalingletscher. This model enables to account for various geometric config- 
urations, interaction between sliding and tension cracking and water flow at 
the bedrock. We could show that both a critical geometrical configuration of 
the glacier tongue and the existence of a distributed drainage network were 
the main causes of this catastrophic break-off. Moreover, the analysis of the 
modeling results diagnose the phenomenon of recoupling of the glacier to its 
bed as a potential new precursory sign announcing the final break-off. This 
model casts a gleam of hope for a better understanding of the ultimate rupture 
of such glacier sliding instabilities. 



1. Introduction 

Gravity- driven instabilities include landslides, mountain collapse, rockfalls, ice 
mass break-off and snow avalanches. They pose a considerable risk to mountain 
communities, real-estate development, tourist activities and hydropower energy 
generation. Gravity-driven instabilities are the most widespread natural hazard on 
Earth. In the US and Europe, they are particularly significant and cause billions of 
dollars and euros damage each year. Prediction of such gravity- driven instabilities 
could be used by policy makers to decide on possible evacuation of the potential 
dangerous area. Unfortunately, accurate prediction of such phenomena remains a 
challenge. 

In this context, studying glacier break-off is of particular interest because a 
glacier consists of a unique natural material (ice) where the interface between ice 
and bedrock is well defined. This relative simplicity of the system allows to particu- 
larly focus our study on the rupture processes leading to the initiation of the insta- 
bility. Two types of glacier instabilities leading to icefalls may be distinguished: (i) 
mechanical instabilities and (ii) sliding instabilities. Mechanical instabilities mainly 
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affect unbalanced cold hanging glaciers (i.e. snow accumulation is mostly compen- 



sated with break-off, Prolong and Funk (2006)). Based on field data combining 



surface displacement measurements and seismic activity before break-off, Faillettaz 



et al (20081 |2011a| ) showed that a prediction of the event was possible. The case 
of sliding instabilities is more problematic. They occur on temperate glaciers (i.e. 
glaciers that can slide on their beds) where subglacial hydrology plays a major role. 
Melt water flowing on the glacier bed may influence the glacier dynamics in two 
ways: it allows (i) a better lubrication of the bed and (ii), if water becomes pres- 
surized, an uplift decouples the glacier from its bed. In the Alps, only two glaciers 
are known to have given rise to such sliding instabilities in the 20th century: The 
glacier du Tour (Mont Blanc, France) in 1949 and the Allalingletscher (Valais, 
Switzerland) in 1965. Up to now, such glacier instabilities are still unpredictable, 
and the causes of such catastrophic break-off remain unclear. 

The present paper is devoted to the study of such instabilities, taking the Allal- 
ingletscher as an example. This glacier is indeed of particular interest because it 
experienced 2 catastrophic break-offs (in 1965 and 2000) and also because a unique 
set of data was collected since the first catastrophic break-off in 1965. To address 
the open questions on the initiation of the instability, we apply a general numerical 
model developped to investigate gravity- driven instabilities in heterogeneous me- 
dia. This model was already successfully applied to a polythermal glacier becoming 



partly temperate at its bedrock, i.e. the Altelsgletscher iFaillettaz et al, 2011b). 



Numerical results showed that the final instability is driven by a progressive cold- 
temperate transition at the ice/bed interface. The present study significantly differs 
from the previous case as here the whole glacier is temperate and can slide on its 
bedrock. As subglacial water was shown to drive the sliding instabilities, the initial 
model was extended to account for subglacial water flow at the bedrock. 

After describing the Allalingletscher and analysing available measurements, we 
apply a general numerical model to this particular gravity-driven instability to as- 
sess the causes of the rupture. Numerical results are discussed and general criterions 
leading to the instability are presented. 

2. Allalingletscher 

2.1. Allalingletscher through the centuries. Allalingletscher is referred to as 
early as 1300 in written documents, in connection with the right of pasturage. The 
glacier is mentioned again and again, through the centuries, mainly in reference 
to floods caused during stages of glacier advance, and also when it interrupted the 
path leading to alpine meadows and pass routes to Italy. The hydrological bassin 
of Allalin (12 km 2 , 81 % glacierized) , located near the head of the valley of Saas, 
consists of a hanging valley, from which the lower part of the glacier tongue extends 
and flows down the side wall of the main valley. During the little age, the terminus 
reached and blocked the river Saaser Vispa leading to the formation of the lake of 
Mattmark. 

The extend of the Allalingletscher has been inferred by Liitschg (see reference 
in Rothlisberger and Kasser (1978)) from references to the size of the lake and to 
the floods which originated in the outburst of the ice-dammed lake. Combining 
such information with direct references to the size of the glacier, he assembled a 
detailed chronicle of the glacier through the centuries, listing large extension of 
the glacier in 1300 and 1589, and on other occasions during the 17th, 18th and 
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Figure 1. Allalingletscher in 1916 (a) and in 2004 (b) 



19th centuries. Maximum known extends occurred in 1822/1823 and 1858. The 
20 th century maximum in 1923 was smaller, although the terminus still crossed the 
main river. After that, an almost uninterrupted retreat lasting 31 years from 1923 
(Fig. [l| ) until 1954 occurred. According to the records published in Glaciological 
reports ( 1881-2009[ ), the retreat amounted to about 230 m from 1923 to 1943, 660m 
from 1943 to 1949, and 90m from 1949 to 1954. The very rapid retreat in the forties 
was due to the decay of the thinned-down snout on the east-facing rock slope of 
the lower part of the Saas valley sidewalk Between 1954 and 1964 the terminus has 
been approximately stationary (Fig. [2|. 

The flat area upstream the south lateral morain was used for the construction 
of an artificial reservoir. Between 1958 and 1967 an earth dam 120m high with 
a volume of 10 7 m 3 was erected immediately upstream the south of the moraine. 
A major part of the material used in the construction of the dam was obtained 
from the lateral moraines, while the flat area between the moraines was used for 
temporary installations, primary in connection with the processing of the earth-fill 
materials. 



2.2. The ice avalanche of 30 August 1965. The Allalin ice avalanche, a landslide- 
type event, is also referenced to as the Mattmark glacier catastrophe. Approxi- 
mately 2.10 6 m 3 of ice broke off at the terminus of Allalingletscher, moved down 
a rock slope of some 27° over the vertical distance of 400 m and continued for 
a further 400 m across the flat bottom of the valley, claiming 88 victims at the 
Mattmark construction site. An overall view of the area shortly after the avalanche 
is given in Fig. [3j 

Glaciological investigations showed that the ice avalanche has occurred during a 
phase of enhanced basal motion as a result of in tensive bed-slip of an even large r 
mass than the one that broke off on August 30 th (Rothlisberger and Kasser 1978). 



Although it seemed that this particular condition was necessary for the lower part 
of the ice mass to slide off, it does not explain why the ice avalanche occurred in 
1965 and not during other active phases either before or after this date. A certain 
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Figure 2. Digital Elevation Model and different extensions of the 
glacier tongue since 1860. 




Figure 3. The Mattmark catastrophe: 30 August 1965 Breaking 
off of 2 mio m 3 of ice 



topography of the bed, combined with an unfortunate mass distribution, is believed 
to have played a major role in the catastrophe. 

In terms of glacier variations, the avalanche results in a sudden retreat of the 
terminus by 400 m in plan view and 220 m in elevation. The avalanche deposit was 
not counted as part of the glacier in this context, although it took more than a year 
for the last remnants of the debris to melt completely. The new terminus consisted 
of a 40 m high ice cliff in the form of a concave crescent-shape arch. 

2.3. Glaciological investigations after the ice avalanche (1965-1978) after 
Rothlisberger and Kasser ( |1978| ). An intensive glaciological study was under- 
taken immediately after the catastrophe in order to safeguard the rescue operation, 
and also the construction work when it was resumed later on. Since the termina- 
tion of the construction of the dam in 1967 a reduced program of measurements 
has been carried out. 

One of the methods of investigation consisted of aerial surveys repeated at short 
intervals of duration between one to two weeks. The first flight was carried out on 
4 September, i.e. 5 days after the fall of the avalanche. 

A detailed monitoring of surface velocities was carried out by means of a theodo- 
lite from the rock ridge of Schwarzbergkopf at the south side of the glacier and stakes 
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drilled in the ice fall. One of the practical difficulties was the installation of stable 
survey markers on this steep and highly crevassed glacier tongue. 

Above the starting zone of the ice avalanche of 1965, a section of medium steep- 
ness remained downstream the flat upper part of the tongue. It consisted of two 
parts, namely an ice-fall with towering seracs separated by deep crevasses followed 
by a detached ice mass which did not fall down. This remaining ice mass expe- 
rienced enhanced motion revealed by several features. Firstly, the ice mass was 
separated from the glacier by a crescent-shaped depression full of shattered ice. 
Secondly, numerous sharp-edged, obviously fresh crevasses cut the inner part of the 
mass involved, which further was more intensely shattered along one edge of the 
glacier. Thirdly, the talus at the foot of the glacier front was pushed over the edge 
of a terrasse, a clear sign that the ice cliff was undergoing rapid advance. 

In the photographs of the construction site with the Allalingletscher in the back- 
ground it no fresh crevasses could be observed on August 5 but, by August 30, clear 
ones already existed. This indicates that the catastrophic event occurred during a 
phase of enhanced motion, called active phase. 

One of the most striking facts illustrated by Fig. [4] is that the timing of the 
beginning of the active phase was different each year, although its end occurred 
at the same time in 1965 and 1966. In 1965 the velocity was already fairly high, 
although still rising, when the survey started. At the beginning of November 1966, 
the velocity record was interrupted when the work on the dam was stopped because 
of the high surface velocity; most markers had disappeared when the observations 
were resumed at the end of the month. In 1967 the survey was discontinued before 
the active phase ceased completely, but the velocity records of the three years 
nevertheless gives a unique picture of the active phase. A common feature in 
the different years is a regular increase and decrease of velocity with time. The 
regularity of the velocity increase is particularly interesting, since a similar pattern 
has been observed on other glaciers prior to the breaking-off of a large ice masses; it 
is noteworthy, however, that the active phases ceased suddenly without formation 
of any major avalanches similar to that of 1965. 

2.4. Glaciological survey since 1978. After the intensive glaciological surveys 
performed by Rothlisberger and Kasser (1978) after the 1965 event up to 1978, 



aerial surveys were performed each year. Complete stereophotogrammetry surveys 
were also performed in 10 years interval to generate a Digital Elevation Model 
of the glacier surface topography with 25 m of mesh grid (Fig. [5|. Moreover, 
stereophotogrammetry surveys were performed every year since 1988 to investigate 
the ice thickness change of the glacier tongue along four different flow lines (see 
Fig. [6|. Mass balance of Allalingletscher was reconstructed with the model GERM 



(Glacier Evolution Runoff Model) developed by Huss et al. (2010). An enhanced 



version of this model developed by \Farinotti et al (2011) was used to determine 



daily runoff for Allalingletscher catchment. 

Thanks to available comprehensive database, it was possible to investigate the 
geometrical evolution of the glacier tongue, which is shown in Fig. (7^a). According 
to this figure, the glacier readvanced rapidly after 1965 and recovered its previous 
geometry 5 year later. The glacier continued to advance up to 1984. Since then, 
the glacier retreated, with an accelerating rate after 1997. 

In August 2000, the glacier configuration was similar to 1965, and for security 
reasons the hazard zone was closed during summer. An ice volume of 1 Million 
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Figure 4. Evolution of the surface velocity of Allalingletscher 
tongue between 1965 and 1967. 
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cubic meters broke of, but did not cause any damage (Fig. [8]). After this event, 
the glacier was stable, and the terminus stayed more or less at the same altitude. 
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Figure 6. Profiles along which glacier thickness was evaluated. 



2.5. Initiation of active phases. Seasonal changes in surface velocity have been 
observed in the years after the catastrophe, and it is now known that Allalin- 
gletscher speeded up regularly every 1-3 years, usually during summer or late au- 
tumn. In most cases no large release of ice occurred. Apparently the active phase 
is necessary but not sufficient to cause breaking off. 

Figs. [7] also emphasizes the years when active phases were observed (red dots). 
At a first glance, they seem to appear when large rate of change of the tongue 
geometry occur, as they are concentrated between 1965 and 1980 (during advance) 
and between 1994 and 2000 (during retreat). However, no clear correlation could 
be found, suggesting another initiation process. The active phase has always begun 
during the melt season, and, perhaps even more significantly, has stopped at the 
beginning of the winter. This seems to indicate that meltwater plays a major 
role in the processes triggering the active phases. However, the initiation of such 
phenomenon does not seem to be correlated with the total amount of water flowing 
each year under the glacier (Fig. (7^c)). 

Ice thickness changes at the tongue was analyzed every year since 1990 by in- 
terpolating the measured thicknesses along the 4 profiles given in Fig. [6j Fig. [9] 
shows the difference in tongue thickness from one year to the next. 

Fig. [7ja) shows that, during this period, active phases occurred in 1994, 1996, 
1997, 1998 and 1999. Interestingly, between 1990 and 1993, the geometrical exten- 
sion of the glacier tongue was stable and the thickness of the tongue experienced 
a small homogeneous thinning. No active phases were detected during this period. 
Then, in 1993, the glacier started a progressive retreat and a thickenning of the 
upper part of the tongue could be observed. This initial "pulse" is due to the pos- 
itive mass balance in 1993 (Fig. (7^b)). In 1994, the tongue became thinner in the 
upper part while a thickening was observed behind the terminus. An active phase 
was detected in 1994. Since then, a "oscillating behavior" was initiated, where ice 
mass is first accumulated in the upper part of the tongue and then transferred to 
its terminus the following year. Such oscillating behavior and active phases seems 
to be correlated. 
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Figure 7. (a) Cumulative length change of Allalingletscher since 
1880. Red dots indicate when active phases occurred, (b) Recon- 
structed mass balances since 1965. (c) Evolution of the mean daily 
runoff of subglacial water flow per each year. 
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2.6. Summary of the Allalingletscher break-off. 

(1) There are two known events: 31.8.1965 and 31.7.2000. 

(2) The glacier tongue is temperate: the glacier is sliding on its bedrock. 

(3) A regular speed up of the glacier tongue during 2-3 weeks has been observed 
almost every year after 1965 between July and October. 

(4) Subglacial hydrology plays a major role on the instability. 

(5) The speed up of the glacier tongue is only a necessary condition for the 
breaking off. 

(6) A critical mass distribution within the glacier tongue is probably a key 
factor for the instability. 

To investigate in more details the causes of this instability we reanalysed this 
event by applying a new numerical model designed for describing natural gravity- 



driven instabilities (Faillettaz et al , 2010). This model allows us to test the different 



hypothesis proposed previously to explore the possible causes of the break-off of this 
glacier. 



3. Numerical modeling 

3.1. Model description. We use a model describing the progressive maturation 
of a mass towards a gravity-driven instability, which combines basal sliding and 
cracking. Our hypothesis is that gravity-driven ruptures in natural heterogeneous 
materials are characterized by a common triggering mechanism resulting from a 
competition between frictional sliding and tension cracking. Heterogeneity of ma- 
terial properties and dynamical interaction of damage and cracks along the sliding 
layer seem to have a significant influence on the global behavior and have to be 
modeled. 
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Figure 9. Difference of thickness of the glacier tongue between 
two successive years from 1989 to 2000. Spatial extension of the 
tongue is also shown. 



This numerical model is based on the discretization of the natural medium in 
terms of blocks and springs forming a two-dimensional network sliding on an in- 
clined plane (Fig. 10 Each block, which can slide, is connected to its four neighbors 
by springs that can fail, depending on the history of displacements and damage. 
We develop physically realistic models describing the frictional sliding of the blocks 
on the supporting surface and the tensile failure of the springs between blocks 
proxying for crack opening. Frictional sliding is modeled with a state-and-velocity 
weakening friction law with threshold. This means that solid friction is not used as 
a parameter but as a process evolving with the concentration of deformation and 
properties of sliding interfaces. Crack formation is modeled with a time-dependent 
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Figure 10. Illustration of the model consisting of spring-blocks 
resting on an inclined plane. The blocks lie on an inclined curved 
surface and gravity is the driving force. Only a small subset of the 
spring-block system is shown here. 



cumulative damage law with thermal activation including stress corrosion. In or- 
der to reproduce cracking and dynamical effects, all equations of motion (including 
inertia) for each block are solved simultaneously. 

The present model improves the multi-block model of \Andersen et al\ ( 1997| ) 



and \Leung and Andersen (1997) in two ways. First, we use a state-and- velocity 



weakening friction law instead of a constant (or just state- or velocity- weakening) 
solid friction coefficient. Second, rather than a static threshold for the spring fail- 
ures, we model the progressive damage accumulation via stress corrosion and other 
thermally activated processes aided by stress. Both improvements make the numer- 
ical simulations significantly longer but present the advantage of embodying rather 
well the known empirical phenomenology of sliding and damage processes. Adding 
the state and velocity-dependent friction law and time-dependent damage processes 
allows us to model rather faithfully the interplay between sliding and cracking be- 



tween blocks and the overall self-organizing of the system of blocks ( Faillettaz et 



o/J[2010|. 



The geometry of the system of blocks interacting via springs and with a basal 



surface is depicted in Figure 10 To sum up, the model includes the following 
characteristics: 

(1) Frictional sliding on the ground, 

(2) Heterogeneity of basal properties, 

(3) Possible tension rupture by accumulation of damage, 

(4) Dynamical interactions between blocks along the sliding layer, 

(5) Geometry and boundary conditions, and 

(6) Interplay between frictional sliding and cracking. 
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We describe the two key processes in the model, the friction and damage laws, 
that are applied to blocks and bonds respectively. 

3.2. Friction law between the discrete blocks and the basal surface (after 



Faillettaz et al. (2011b). Guided by the rate-state friction character of rock 
~( Dieterich\\1994§ and ice \Fortt and SchulsorA 120091 \Lishman et all |2011|) , a rate- 



and state-dependent friction law seems also to be adequate to describe ice-bedrock 
friction. The version of the rate/state- variable constitutive law, currently most 
accepted as being in reasonable agreement with experimental data on solid friction, 



is known as the Dieterich-Ruina law Dieterich (1994): 



(3-1) 



/i(i,0) = /i + Aln^- 



Bin- 



where the state variable is usually interpreted as the proportional to the surface 
of contact between asperities of the two surfaces. The constant /iq is the friction 
coefficient for a sliding velocity So and a state variable 0q. A and B are coefficients 
depending on material properties. 

The time evolution of the state variable 6 is described by 



(3.2) 



dO 
dt 



= 1 



0S_ 



where D c is a characteristic slip distance, usually interpreted as the typical size 



of asperities. The friction law (3.1) with (3.2) accounts for the fundamental prop- 



erties of a broad range of surfaces in contact, namely that they strengthen (age) 
logarithmically when aging at rest, and tend to weaken (rejuvenate) when sliding. 

The primary parameter that determines stability, A — B, is a material property. 
Physically, A — B indicates the sensitivity of the friction coefficient to velocity 
change: a negative value indicates velocity- weakening leading to an unstable slip; 
a positive value indicates velocity-strengthening leading to a stable slip. Faillettaz 
et al. (2010) showed that the case A = B is of special interest because it retrieves 



the main qualitative features of the two classes, and also because, empirically, A is 
very close to B. Assuming A = B is therefore reasonable and ensures more robust 
results. In Appendix A we present how the equation for the critical time tf was 
obtained, tf indicates the transition from a subcritical sliding to the dynamical 
inertial sliding. It is given in this case (Faillettaz et al. 2010) by 



(3.3) 



exp(™) - 1 



for fi > /io; Note that tf — » oo for \i < /xq, where tf is the time when the block starts 
sliding, fiQ is a constant friction coefficient, A is a constant parameter depending on 
material properties and 0q is the state parameter at steady state. The parameter 
\± is evaluated for each block with the definition of the solid friction ji = ^ where 
T is the tangential force determined from the position of its connected neighbors 
and N is the normal component of the weight of the block. The state parameter #o 
in Eq. |3.3| is given (see Appendix A) by 



(3.4) 



0o = 



Dc 
So 
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where So is generally interpreted as the initial low velocity of a sliding mass before 
it starts to accelerate towards its dynamical instability and D c can be interpreted 
as a characteristic slip distance over which different asperities come in contact. 

Three parameters have to be correctly evaluated to model the frictional processes 
within the glacier: /i , A and 0q. 

An increasing friction coefficient with increasing sliding velocity (velocity-strengthening) 
at low velocities has been observed for polished ice-on-ice ( Kennedy et al. t 20001 



Montagnat and Schulson, 2003), ice-on-granite [Barnes et al. 1971) and ice-on-ice 



along Coulombic fault \Fortt and Schulson\ |2009| ). A decreasing friction coefficient 
with increasing sliding velocity (velocity- weakening) at high velocities has been ob- 
served as well ( Barnes et al.\ |1971| \Kennedy et al\ |2000). The transition between 



these two regimes occurs at a velocity of 10 ms (Fortt and Schulson, 2009). 



These two different behaviors are generally explained by two physical mechanisms 
depending on the sliding velocity regimes: 



(i) 



The first mechanism is the water lubrification mechanism (produced by 
frictional heat at the sliding surface) working at sliding velocities above 
roughly 0.01 m s _1 . The water lubrification mechanism is characterized by 
the low viscous resistance of a water film produced by frictional heat at the 



(ii) 



sliding interface ( Maeno and Arakawa 



et al 2000 Montagnat and Schulson 



2003). 



2004 Barnes et al , 1971 Kennedy 



The second mechanism is the adhesion and plastic deformation of ice at 
the friction interfa ce, which is present at velocities lower than roughly 0.01 
m s _1 ^Kennedy et 



Arakawa 2004) 



2000 Montagnat and Schulson, 2003 



Maeno and 



By analogy with rock physics, a rate- and state-dependent friction model for 



fresh and sea ice was proposed recently ( Fortt and Schulson 2009 



Lishman et al. 



2011). Such a model enables to explain the change in the friction coefficient when 



varying sliding velocity along Coulombic shear fault or slip history. The experiments 
performed by Fortt and Schulson] ( |2009| at -10°C on sliding along Coulombic shear 
faults in ice suggests that the ice-ice friction coefficient is velocity- dependent and 
varies between 0.6 to 1.4 depending on the applied sliding velocity along the fault. 
As these last results are the closest to our case, we used a mean value of ii l ^ e ~ lce 
equal to 0.8 in our calculation. 

Once the block slides, the dynamics is controlled by a kinetic friction coefficient, 
which is in general smaller than the static coefficient fi l ^ e ~ lce . For low temperatures, 
Kennedy et al. (2000) and Fortt and Schulson (2009) found a relatively constant 
friction coefficient around 0.6 for sliding velocities below 10 _4 ms _1 , and rapidly 
decreasing values down to 0.1 at velocities above 10 _4 ms _1 . Because velocities 
above 10 _4 ms _1 are not expected in our model describing the nucleation phase to 
the catastrophic rupture, we thus assume a kinetic friction coefficient Hd = 0.6. 

Next, we turn to the evaluation of the coefficient A (and B since we have made the 
simplifying assumption that B = A) in the rate- and state-dependent friction law 



(Eq. 3.1). Laboratory experiments suggest that A is smaller than fiQ by typically 



one and sometimes up to two orders of magnitude for rock \ Ohmura and Kawamura\ 
|2007| \Scholz 2002, 1998), ranging from 0.01 to 0.2. A recent experiment on saline 
ice suggests however that A could be significantly higher for ice (i.e. A « 0.3) 



(Lishman et al. 2011 ). As we do not have access to strong experimental constraints 



and since the friction law should be valid for the ice-ice and the ice-bedrock friction, 
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we choose A « 0.1, corresponding to one-tenth of the static friction coefficient, 



which is the upper limit given by Ohmura and Kawamura (2007) for rock. 



The last parameter to be determined is 0q m Eq. |3.4| In the case of a glacier, the 
sliding velocity is typically of the order of centimeters per day. Therefore a roughly 
correct estimation is So ~ 1 cm.d -1 . D c can be interpreted as a characteristic slip 
distance over which different asperities come in contact. It is difficult to evaluate 
this value. The recent seismological literature reports D c to lie in the range of tens 



2003 Zhang et al. , 2003). 



1 m. Finally, inserting this value in equation (3.4), we 



of centimeters to meters for earthquakes ( Mikumo et 
We arbitrary choose D, 
obtain 6q = 100 days. 

To account for the heterogeneity and roughness of the sliding surface, the state 
variable 6i is reset to a new random value after the dynamical sliding stops. This 
random value should not be choosen too low in order to prevent a block which just 
stopped sliding from switching immediately in a new dynamical phase (i.e. tf = 0). 
Thus, we assign 0i = v.Qq with v uniformly distributed between 0.5 and 1.5. 

Finally, the possible water lubrication of the bedrock leading to a progressive 
decrease of the ice-rock frictional resistance is modeled using a progressive decrease 
of the friction coefficient /iq at the interface between the simulated ice block and 



the basal surface (see section 3.5). 



3.3. Creep law. As explained in section [3TT[ bonds are modeled as linear springs 
in parallel with an Eyring dashpot. The springs transmit the forces associated 
with the relative displacements of the blocks. The spring stiffness has also to be 
evaluated in order to reflect the elastic property of the bulk ice mass. In continuous 
elasticity, Hooke's law of elasticity relates stress a and strain e via Young's modulus 
E: a = Ee. This leads to the expression: (Jbuik = ^ice^r- This stress is applied 
to a side surface of the block S = L x H (where H corresponds to the height 
and L the length of the surface), leading to an equivalent force in the bulk equal to 
Fbulk = Eice^t S - A nnear spring is subjected to forces given by F bond = K bond SL. 
In order to find an equivalent behavior, these two forces have to be of the same 
order, leading to a spring stiffness given by K bon( i = Ei ce H. Usually, values for 
Ei ce are reported to be 9 GPa (Petrovic t 2003 Petrenko and Whitworth 1999). 



However, there is a disagreement of an order of magnitude between measurements 
of E in laboratory (9 GPa) and from field observation lGPa) as argued by 
Vaughan\^mE ). 

Depending on the applied stress and the time scale of interest, ice has either 
a linear viscous or a non-linear viscous behavior. In glaciers, ice creep is usually 



described with a non-linear viscous rheology called Glen's flow law (see Hutter 



(1983) and references therein). This law relates, in steady-state conditions, strain 
rate and stress in the secondary creep regime. It is thus not possible with this law 



to describe the tertiary creep and the rupture of a bond. However, Hutter (1983) 



states that the so-called Prandtl-Eyring flow model shows a compatible behavior 
with Glen's flow law at low stresses. Non-linear viscous behavior is introduced 
in our model with an Eyring dashpot in parallel with a linear spring of stiffness 



E (which is analogous to the Prandtl-Eyring flow model), following Nechad et al. 
( 2005 ) . Its deformation e is governed by the Eyring dashpot dynamics and reads 



(3.5) 



de 

— = i^sinh(/35 dashp ot), 
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where the stress Sdashpot m tne dashpot is given by 

s 



(3.6) 



"^dashpot 



1 - P(e) 



Ee . 



Here, s is the total stress applied to the bond and P(e) is the damage accumulated 
within the bond during its history leading to a cumulative deformation e. P(e) can 
be equivalently interpreted as the fraction of representative elements within the 
bond which have broken, so that the applied stress s is supported by the fraction 
1 — P(e) of unbroken elements. Following \Nechad et al. (2005), we postulate the 
following dependence of the damage P(e) on the deformation e: 

(3.7) P(e) = 1 - 



eoi 



e 2 



where eoi, eo2 and £ are constitutive properties of the bond material. 

( [2010D ended up 



Finally, by combining the previous equations, Faillettaz et al. 



with a creep model that computes the critical time (i.e. failure of the bond) as a 
function of the stress experienced by the bond s given by: 



(3.8) 



where 



(3.9) 



and 



(3.10) 



K 



exp(— 7 s) 



/?4 



-oi 



if s > 
if s < 



E 



(?) ( 



eo2 



Creep properties are defined by the parameters K, (3, eoi, eo2 and £, that we need 
to fix for our simulations. 

We need to find the most appropriate parameters to describe creep behavior of 
ice. The ice of natural glaciers has a complex polycrystalline structure composed 
of crystals of different sizes. As Equation 3.7 shows, a fraction 1 — (eoi/eo2)^ of 
all present representative elements undergo abrupt failure immediately after the 
stress is applied. Since ice does not show significant damage immediately after 
being loaded, a reasonable assumption is eoi = eo2- The relative heterogeneity of 
the material is introduced with the parameter £. Despite its complex structure, ice 
is a fairly homogeneous material compared to fiber matrix composite. The more 
homogenous a material, the greater £. In the following, we set £ = 10 (which 
means that ice is not very heterogeneous, for example, Nechad et al. ( 2005| ) 
£ ^ 2 for fiber matrix composite). Moreover, as eoi 



used 



eo2, £ has just an influence 
on the parameter s*, i.e. the critical stress above which damage starts. 

The other parameters describing the deformation of the Eyring dashpot under 
an applied stress are j3 for the non-linear term and K for the linear one. It is 



difficult to infer such parameters for ice. Nechad et al.\ ([2005]) used for fiber matrix 

1 (F~ 



composite f3 = 50 10" 9 Pa _i and K = lir s 
arbitrary choose f3 = 10~ 7 Pa -1 and K = 10 
ice equal to 1 MPa, we obtain from Eq. 3.5 ^ = 



the behavior of polycrystalline ice ( Schulson and Duval 2009) 



As ice is significantly weaker, we 
3 s _1 . With a tensile strength of 
: 10 _4 s _1 , which is coherent with 
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Slope distribution under blocks 




Slope (degree) 

Figure 1 1 . Distribution of the slope of the bedrock at the position 
of the blocks. 



3.4. Geometric parameters. We have to consider the geometric input parame- 
ters for modeling Allalingletscher. The glacier is discretized into a system of regu- 
lar cubic blocks. As we saw in section |2.4[ spatial extension of the glacier tongue 
evolved a lot during the last century. Digital elevation models (DEM) of the tongue 
for different years are available (see Fig. |5|. In this study, we perform simulations 
for two extreme cases, one where the extension of the glacier tongue was maximum 
i.e. in 1982 and the other where the glacier broke-off i.e. in 2000. Since the glacier 
has completely retreated from the steep part, a detailed bedrock topography of the 
former glacier tongue could be obtained by a photogrammetric processing of recent 
aerial photographs (2008). 

In order to obtain a realistic description of the damage and fragmentation process 
that may develop in the ice mass, we need a sufficiently large number of blocks. 

As a compromise between reasonable sampling and numerical speed, we use a 
model composed of 20 m x 20 m blocks. 

The spatial extensions of the glacier as well as the block heights in 1999 and 
1982 is given by the respective DEM. The slope <j) of the bedrock ranges from 0° to 



45° (see Fig 11). 



3.5. Extension of Faillettaz et al. (2010) model: accounting the effect of 
subglacial water flow. In our initial model, the effect of subglacial water flow on 
basal motion was not accounted. Basal processes play a key role in the instability 



(see|2.5). Melt water flow at the glacier bed influences the glacier dynamics ( Be, 
2008) in two ways: (i) it lubricates the bed and (ii), if water becomes pressurized, 



the glacier can be partially decoupled from its bed. The lubrication of the bedrock 
leading to a progressive decrease of the ice-bedrock frictional resistance is modeled 



by decreasing the friction coefficient /io ( m Eq. 3.3) {Faillettaz et al, 2011b). The 
following assumptions are made: 



(1) The friction coefficient depends on the subglacial water pressure {Schweizer 
and Iken t 1992 Jay-Allemand et al. 2011). 
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(2) The subglacial water pressure and discharge variations at the daily time 



scale are similar (Boulton et al. 2007). 



(3) As a consequence of the two aformationned assumptions, the decrease of 
the friction coefficient can be assumed to be proportional to the discharge 
of subglacial water flow. 

(4) The ice thickness does not vary drastically along the glacier tongue. Hence 
the gradient of the hydraulic potential and therefore the flow path of water 



is directly given by the topography of the bedrock (Flowers and Clarke 



1999) 



At a first step, the amount of basal water flowing under each block has to be 
evaluate. We used a toolbox developped by Schwanghart and Kuhn (2010|) to assess 



the subglacial drainage network with a DEM of the glacier bed. Most commonly 
DEMs are represented as rectangular grids where an elevation value is assigned to 
each cell. Two steps are needed to evaluate the basal water flow under each cell of 
the DEM. First, the flow direction is evaluated via a transfer matrix M. Elements 
in M contain the relative amount of discharge Mij transferred from one cell i to 
a maximum number of eight downward neighbors with the index j. Thereby the 
transfer ratios are proportional to the downward slope to the respective neighbor. 
Second, by counting the number of cells draining in each grid cell, it calculates 
the upstream contributing drainage area. Finally, this toolbox gives the relative 
contribution from each cell of the DEM to the global runoff at the outlet. The 
amount of water flowing through a particular cell of the DEM is then obtained by 
multiplying the total runoff at the outlet by its relative contribution. Results are 



shown in Fig. 12 



In a second step, the friction coefficient is adapted according to the basal water 
discharge under each block with the following relationship: 

(3.11) Voihj) = Mi - Z ev (i,j) - Q-Cp, 

where fi is the friction coefficient of a given block (i,j), /ii a constant friction 
coefficient, Z ev (i,j) is the relative contribution of the block, Q is the daily runoff at 
the glacier terminus and c p is a parameter describing the relative influence of the 
basal water on the friction coefficient. In this way, the friction coefficient of a block 
is assumed to decrease when the local subglacial water discharge (and therefore the 
local subglacial water pressure) is increasing. 

3.5.1. New algorithm accounting for subglacial basal water flow. The different steps 
describing how the instability is modeled are plotted in Fig. [l3]and developped in 
Annexe 2. As explained previously, two phases have to be distinguished: 

(i) A quasi-static (quiescent) phase corresponding to the nucleation of block 
sliding and bond rupture. 

(ii) A dynamical (active) phase corresponding to the sliding phase of the blocks 
and the failure of bonds. 

The new algorithm takes into account the subglacial water flow. At each time 
step of the run, the friction coefficient under each block is modified according to the 



local subglacial discharge. It is evaluated with Eq. 3.11 based on the runoff at the 



outlet given by the model of Farinotti et al. (2011[) and the matrix Z ( 



from the DEM of the bedrock topography. 
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(a) 



(b) 



Figure 12. Relative contribution of each block to the total sub- 
glacial water discharge (Z ev ) for 1999 (a) and 1982 (b). Red zones 
indicate where discharge of subglacial water is large. 

4. Numerical results 

The aim of the numerical simulations is to better understand the development 
of a glacier sliding instability, considering Allalingletscher as an example. In par- 
ticular, we intend to provide answers to the following questions: 

• Why did the glacier only broke off in 1965 and 2000? 

• Does the glacier geometry play a role in the occurrence of the break-off? 

• What is the influence of subglacial water on the glacier dynamic and on the 
break-off? 

• Does the subglacial drainage network efficiency play a role in the onset of 
the instability? 

• Is it possible to find precursory signs of the break-off? 

• Would it be possible to predict the break-off? 
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Figure 13. Schematic flowchart of this modified spring-block 

model accounting for subglacial water flow. 



4.1. Tuning the creep parameters. The first step of the numerical simulation is 
to find plausible creep parameters able to reproduce the observed glacier behavior. 
Creep parameters turn out to be crucial. As explained in section |3.3[ the creep 



parameters are f3 and C. Faillettaz et al (2010) showed that three different regimes 
can be expected from the model, namely a fragmentation regime, slab regime or 
stick and slip regime. As showed in section [Z4} the glacier broke off in 2000 as a 
slab avalanche. So we have to find a parameter set able to reproduce the observed 
behavior. 

In this part, we will use the 1999 surface topography for our simulation. At each 
time step, we evaluate, for each block, the basal slope (see Fig. [2|. We performed 
different runs with the 1999 glacier geometry by varying creep parameters. Results 
are shown in Fig. [l4j 

In the fragmentation regime, creep occurs before the nucleation of sliding could 
develop, so that bonds break and the bottom part undergoes a fragmentation pro- 
cess with the creation of a heterogeneous number of sliding blocks (Fig. [l4^b)). 
In the stick-slip regime, the time needed for damage initiation is so great that the 
whole glacier tongue undergoes a series of internal stick slip events, associated with 
an initial slow average downward motion. After a sufficient time, all blocks slide 
(Fig. |l4^d)). In the slab regime, in which neither damage nor frictional sliding 
dominates, the occurrence of a macrocrack propagating roughly along the location 
of the largest curvature associated with the change of slope from the stable fric- 
tional regime in the upper part to the unstable frictional state in the lower part 

(Fig. [njc)). 

In the following, we used the parameter set that reproduces the last regime. 



4.2. Influence of the glacier geometry on the instability. To investigate the 
contribution of the glacier geometry on the onset of the instability, we performed 
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Figure 14. Snapshots illustating the three possible regimes: frag- 
mentation (b), slab (c) and stick-slip(d) 

Table 1. Numerical values of parameters used in our simulations. 







Friction parameters 


A 


#0 


Mo #o 




d 


m.d~ l 


0.1 


100 


0.8 IO" 3 






Creep parameters 


E 


p 


C~l/K £ e i eo2 


Pa 


Pa" 1 


s - 


10 9 


io- 7 


10 3 10 0.003 0.003 






Runoff parameters 


Mi 






0.8 


6 





numerical simulation on two different glacier geometries in year 2000 (when a break- 
off occurred) and year 1982 (when active phase could be evidenced without break- 
off). We initialized our simulation as described in the previous section. Table [I] 
summarizes the parameters used in our simulations. 

Fig. 15 shows the results of simulations for both geometries of 2000 and 1982 
with the same creep and friction parameters. As the terminal part of the glacier 
lies on a average slope of 30°, we imposed the friction coefficient to be a bit lower 
than its mean slope, i.e. /iq = 0.5. 

With this set of parameters, the initiation of the instability could be reproduced. 
At the beginning of both simulations, blocks situated in the steepest part start to 
slide, inducing a zone with high tension at the slope changes. Then, blocks located 
downstream start to slide resulting in the apparition of a compression zone directly 
downstream of this process zone (Fig. [l5^b)). This compression zone grows as 
the instability develops. In the case of the 1982 geometry, the compression zone 
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Figure 15. Snapshot describing the rupture progression for the 
glacier geometry in the years 2000 and 1982, with a constant fric- 
tion coefficient. 



continues to propagate downwards and reaches the terminus where the glacier is less 
steep. The whole glacier is then supported by the lower part of the tongue, which 
stabilizes the glacier and stops the development of the instability (Fig. [l5^e2)). On 
the contrary, in the case of the 2000 geometry, the compression zone reaches the 
terminus where the slope is still large. Nothing supports the advance of the glacier, 



resulting in a global break-off of the tongue (Fig. 15 'el)). 
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Figure 16. (a) Blue blocks indicate where friction coefficient was 
changed, (b) Snapshot at the end of the run. 



The geometrical configuration of the glacier tongue plays therefore a major role 
in the development of the instability. In the case of Allalingletscher, the glacier 
tongue where the terminus is located beyond the position of 1965 is stabilized as 
the slope near the terminus is lower. 

4.3. Role of subglacial water: Channelized vs. distributed subglacial 
drainage network. 

4.3.1. Channelized subglacial drainage network. We pointed out that the subglacial 
water and drainage network play a role in the initiation and the development of 
the instability. The presence of subglacial water could be modeled by decreasing 



the friction coefficient as a function of discharge (see 3.5). 

To investigate the relative contribution of the drainage network at the onset of 
the instability, we performed different type of simulations. We have already shown 
in section 



4.1 



that a global decrease of the friction (i.e. modeling a distributed 
subglacial drainage network) could qualitatively reproduce the break-off event in 
the year 2000. 

To study the effect of a channelized drainage network, we performed a simulation 
where the friction coefficient was decreased locally, in a band roughly situated at 
the center of the glacier (Fig [16^l). The band width was set to 3 blocks, i.e. 60 m. 
This allows to reproduce the effect of a flow path of subglacial water in a restricted 



area. Results are shown in Fig. 16 After the development of the crown crevasse, 
the instability develops, leading to a localized break-off with a triangular shape. 
Moreover, to obtain such a rupture, the friction coefficient had to be set to a very 
low value (/io = 0.2), indicating a stronger stability in the case of a channelized 
network than with a distributed one. A channelized network could not reproduce, 
even qualitatively, the break-off event of the year 2000. 

4.3.2. Distributed subglacial drainage network. We performed two runs with the 
same creep and friction parameters but for two different years (1982 and 2000). 
For these 2 years, modeled daily runoff at the outlet of the glacier are available 
(red lines in Fig. 



according to Eq. 



17). The model allows to redistribute the water under each block 



3.11 (Fig. 12). 



Fig. 18 shows the evolution of the instability for the years 2000 and 1982. In 
2000, after the opening of the crown crevasse and the formation of a compressive 
zone at the glacier tongue (Fig. [l7^a)), the glacier enters in a quiescent regime. 
Even the increase in runoff in May does not affect the dynamic of the glacier 
tongue (Fig. [l7^a)). Suddenly, after 165 days of simulation (in mid- June), the 
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Time (a) 
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Figure 17. Runoff, number of sliding blocks and surviving bonds 
as a function of time for the situation in 2000 (a) and 1982 (b) 



tongue starts to be unstable, resulting in a fragmentation process that propagates 
over the whole tongue in a few days. We stopped the simulation when the glacier 
tongue was completely disaggregate. In 1982, the situation is different as the tongue 
extension is larger. The compression zone propagates downstream until the lowest 
glacier tongue, where the bed slope decreases. An increase in the number of sliding 
blocks can be noticed between mid July and beginning of September (Fig. [Pf^b)), 
indicating an ongoing active phase. Interestingly, this active phase is not directly 
correlated with the runoff magnitude, as it drops in August in the course of the 
ongoing active phase. It also appears that once triggered, the active phase needs 
some time to stop. 
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Figure 18. Snapshots describing the rupture progression and slid- 
ing instability in the block lattice for the years 2000 and 1982. 
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5. Discussion 

With the results obtained with our model, we could highlight the mechanisms 
leading to the sliding instability. The main result is that the active phase is triggered 
by an increase of subglacial water flow, its duration depending on the extension of 
the glacier tongue. 



5.1. Existence of active phases and prediction. An ongoing active phase does 
not necessarily lead to a catastrophic break-off. Some further conditions are re- 
quired to trigger the large scale rupture: 

(i) A critical geometric configuration of the glacier tongue, 

(ii) A distributed drainage network. 

This explains why such instabilities are so rare. 

Our model is able to explain and predict the occurrence of the active phases 
observed on Allalingletscher, but not the final catastrophic break-off. For the situ- 
ation in year 2000, the model could reproduce with a good timing the active phase 
appearing before the final rupture. In the following, as the glacier is completely 
disaggregated, the simulation could not be pursued and the occurrence of the break- 
off event could not be predicted. Moreover, the rapid motion of the disagregated 
glacier likely affected the efficiency of the subglacial drainage network. This process 
could lead to a positive feedback: an enhanced basal motion affects the efficiency of 
subglacial drainage network leading to increase subglacial water pressure and basal 
motion further. This positive feedback induced by the onset of the active phase is 
not modeled. 



5.2. Impact of the runoff evolution on the break-off. Our results show for 
both years 2000 and 1982 that fragmentation phases (drop of surviving bonds high- 
lighted by gray bands in Fig. 17) coincide with periods where the daily runoff is 
decreasing. On the contrary, when the runoff increases, the bond breaking process 
is suddenly stopped. Moreover, the glacier motion reacts to runoff variations with 
a delay depending on its size. Fig. 17 indicates a clear correlation between runoff 
magnitude and the number of sliding blocks with a larger delay in 1982 than 2000. 

These results provide new insights on the maturation process of the instability. 
They suggest the following sequence of ongoing processes: 

(1) The onset of the active phases is induced by an increase of subglacial water 
flow with a delay depending on its size. 

(2) The initiation of the fracturation process starts in a period of decreasing 
runoff, i.e. during a rapid recoupling phase of the glacier onto its bedrock. 
During this phase, an intensive fracturation of the glacier tongue is initiated. 

(3) A catastrophic break-off requires the combination of two opposite phenom- 
ena: first, the glacier needs to be in an active phase with a strong enhanced 
basal motion, and second, this active phase must be stopped abruptly with 
a rapid recoupling of the glacier to its bed. 

The likelihood of this process chain can be verified with observations during the 
years 2000 (Fig. [l7^a)) and 1965 (Fig. 19) when the glacier broke off. It appears 
that, for both years, one month before the rupture, the runoff drastically increases, 
which could led to an active phase. During the 6 days prior to the break-off, the 
runoff dropped from 13 to 5 m 3 s _1 in 2000 and from 14 to 5 m 3 s _1 in 1965. 
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Figure 19. Modeled runoff regime in the year 1965 (The red 
star represents the break-off). 



It appears that the process chain leading to the final break-off in 1965 and 2000 
inferred by our modeling results is confirmed by these observations. 

With these new insights explaining the sliding instabilities for Allalingletscher, 
we can analyse other situations. 

5.3. Analysis of observed sliding instabilities of other glaciers. 

5.3.1. Break-off of Le glacier du Tour (France). The terminal tongue of the glacier 
du Tour (Mont Blanc, France) broke off on August 14th 1949. The estimation of 
the volume of the break-off ranges from 0.5 to 2 millions cubic meters (Glaister, 
1951). At that time, the glacier geometry was similar to Allalingletscher. The 
upper part of the glacier was lying on a flat zone, whereas its tongue ended in 
a very steep terrain (about 40 degrees) that is not able to support the glacier 
tongue. The geometrical configuration of the glacier tongue was also in a critical 
configuration. Moreover, different subglacial streams emerging from the glacier 
could be distinguished, indicating a distributed subglacial drainage network (see 
Fig. 20 ). The two main criteria for a break-off in the case of a sliding instability were 



likely met. Meteorological conditions prevailing in this area were analyzed (Fig. 21 ) 
to infer a runoff history prior to the break-off event. This temperature time series 
indicates that in the month preceding the break-off, temperatures were high (about 
10 °C at 2472 m). Five days before the break-off, the temperatures dropped from 
14 to °C. Subglacial runoff magnitude should have drastically decreased during 
this period. This confirms that the recoupling of the glacier on its bed played a 
decisive role in the break-off event. 

5.3.2. Triftgletscher (Bernese Alps, Switzerland). The Triftgletscher is situated in 
the Bernese Alps in Switzerland. The glacier experienced a rapid retreat during 
the last 20 years. After the year 2000, the tongue ended in a steep terrain (about 
35 degree) and its stability was questionable (Dalban et al. 2011). In the following 
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Figure 21. Temperature and precipitation in summer 1949 at 
the Col du Grand Saint Bernard (2472 m a.s.l.), located about 
20 km South West from the Glacier du Tour. The break-off is 
indicated with a red star. 



summer surface velocities increase from 1 to 4 md _1 , similar to a mini-active phase. 
However, up to now, the glacier did not break-off. 

All the basic conditions for the glacier to break-off were fullfilled, except the 
existence of a distributed subglacial drainage network. This might be the reason 
for the absence of the initiation of active phase here. 

5.3.3. Giesengletcher (Bernese Alps, Switzerland). In a more general context, cli- 
mate change may affect the stability of sliding glaciers. As a general glacier retreat 
in the Alps is observed, the geometry of some glaciers might potentially evolve 
toward a critical situation. As an example, Giesengletscher (Bernese Alps, Switzer- 
land) would be a suitable candidate to a future catastrophic break-off. The termi- 
nus of Giesengletscher is located at about 2500 m a.s.l. in the Bernese Alps (Fig. 



22). In 2008, a crevasse spanning the whole glacier was observed on the glacier 



tongue indicating an ongoing active phase on the steepest part (about 35 degree, 



Fig. 22 'a)). The situation is nevertheless not critical yet, as the glacier terminus is 



resting on a moderate slope, which stabilizes the glacier tongue (green zone in Fig. 
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Figure 22. (a) Situation of Giesengletscher in 2008 (highlighted 
in red), (b) Situation in 2011. Green zone indicates the stable part 
of the tongue, in red the possible unstable zone. 

[22^b)). Moreover, its bedrock is likely smooth, so that the drainage network is ex- 
pected to be distributed. Except for the geometrical configuration of the tongue, all 
conditions for the glacier to break-off are fulfilled, indicating that the disappeance 
of the supporting glacier terminus could lead to a critical situation, as the tongue 
geometry would become critical. 

6. Conclusions 

Glacier sliding instabilities may occur on temperate glacier tongues. Such in- 
stabilities are strongly affected by the subglacial hydrology: infiltrated melt water 
may indeed cause (i) a lubrication of the bed and (ii) a decrease of the effective 
pressure at the glacier bed and consequently a decrease of basal friction. Available 
data from Allalingletscher (Valais, Switzerland) indicate that the glacier tongue 
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experienced an active phase during 2-3 weeks in summer in most years between 
1965-2000 with strongly enhanced basal motion. This glacier broke-off twice in 
1965 and 2000. In order to scrutinize in more detail the processes governing the 
sliding instabilities, a numerical model developed to investigate gravitational in- 
stabilities in heterogeneous media was applied to Allalingletscher. This modified 
spring-block model enables to account for various geometrical configurations of the 
glacier and also for interactions between basal sliding and tension cracking. The 
impact of subglacial water flow on basal motion was included in the model. 

Our results indicate that basically three conditions have to be fulfilled for a break- 
off to occur in the case of a sliding instability: (i) a critical geometrical configuration 
of the glacier tongue is needed, where the glacier terminus rests on a steep slope, 
(ii) the glacier has to experience an "active phase" during which basal motion 
drastically increases, (iii) the subglacial drainage network has to be distributed. 
Moreover, our modeling results point out that the break-off always occurs after a 
sudden phase of recoupling of the glacier on its bedrock, characterized by a drop 
of subglacial water flow. This new precursory sign was observed for both 1965 and 
2000 break-off event. We also showed that all these conditions were fulfilled in 1949 
before the break-off of the glacier du Tour, and moreover, we demonstrated that 
the same drop of subglacial water runoff occurred five days before the catastrophic 
break-off, indicating that this new precursory sign might be decisive. In a more 
general context, climate change may affect the stability of steep glacier tongues, as 
they retreat and may evolve toward a critical geometry. This model casts a gleam 
of hope for a better understanding of the ultimate rupture of such glacier sliding 
instabilities. 

Initiation of sliding for a single block 
This appendix complements Section |3.2| by providing details of the calculation 



of the critical time at which unstable sliding of a given block occurs. Section 3.2 
describes the sub-critical sliding process of a given block interacting via state-and- 
velocity solid friction with its inclined basal surface. When the sub-critical sliding 
velocity dS/dt diverges (we refer to the time when this occurs as the "critical time" 
tf for the frictional sliding instability), this indicates a change towards a dynamical 
sliding regime where inertia (the block mass and its acceleration in the Newton's 
law) has to be taken into account. 

Let us calculate explicitly how the critical time tf is obtained and define its de- 
pendence on the parameters and boundary conditions. Let us call T = || -Fbond — 
^weight %\\ (o r N = ^V weig ht) the total shear (or normal) force exerted on a given 
block, where Fb on d is the force exerted by a neighboring spring bond, and Af we i g ht 
and T weig ht are the normal and tangential forces due to the weight of the block. We 
then have 

T 

(.1) /i= — =tan</>, 

where <j) is the angle of the basal surface supporting the block. Therefore, 

8 

(.2) /i + A In t- + B In — = tan <j> . 



As explained in Section [3^2} A — B is usually very small for natural material, of the 
order of A — B « ±0.02. For the sake of simplicity, we assume A = B. As recalled 
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in Section |3.2[ this choice is not restrictive as it recovers the two important regimes 



( Helmstetter et a/.|| 2 QQ4| ). This leads to 

5 tan0 - fi 

whose solution is 

( a\ k a (k a \ ftan0-/i O 

(.4) 5.6 = (S 60) exp 



Combining Eqs. |3.2| and[^4| we obtain 

• OS S /^tan0-/i O 

(.5) fl = i_ =i_ exp' 



After integration it reads 



(j 6» ) exp(^=^) 

^C 1 d: z >> 



(.6) 61 
and using Equation (£4j), we obtain 

(S ) exp ( 



tan 0— /io 

(.7) 5 = 4 1 



^0 + (1 - 4^ ex p ( tan ^~ w ) )t 

This expression exhibits the usual regimes: a finite time singularity is obtained for 
(<5o @o ) exp ^ tan< ^~w> j > D c . In this case, Expression can be re- written as 

D c S O exp (iSL^ 



(.8) 5 = 

with 

(•9) */ 



D c - 5 O exp j 



D c - So O exp (iMte ) ' 

We can simplify Expression (^8| by using the condition that, for \i = tan</> = /io, 
we should have tf 00. But, for /i = /io, exp ) = 1 anc [ thus, for the 

condition — » 00 to hold, we need 

MO, 

The final expression for the critical time tf signaling the transition from a subcritical 
sliding to the dynamical inertial sliding is, for /i > /io, 

(•11) t f 



exp(^) - 1 ' 



while tf — » 00 for \i < /iq. Note that the dependence on 5q nas disappeared due to 
the relation po] ). 

To summarize, a given configuration of blocks and spring tensions determines the 
value of T = || J^Fbond — ^weight ^| I an d A/" = A" we i g ht and therefore of /i via 
Knowing /i and given the other material parameters #o,Mo anc i A we determine 
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the time tf for the transition to the dynamical regime for that block via Equation 



.11). 



General Algorithm 

The simulation of the frictional process for each given block proceeds as follows: 

(1) A given configuration of blocks and spring tensions determines the value of 
^ = II Z^bond — ^weight^H and TV = Alight for each block, and therefore 
their solid friction coefficient /i corresponds to the ratio T/N. 

(2) At a given time, knowing the runoff at the outlet, /iq is evaluated for a 



given block with Eq. 3.11 



(3) Knowing \i for a given block and with the other material parameters Oo, Mo 
and A for that block, the time tf for the transition to the dynamical sliding 



regime is calculated with expression (3.3). Eq. 3.3 for tf gives the waiting 
time until the next block starts to slide. 

(4) When the block undergoes a transition into the dynamical sliding regime 
at time tf, its subsequent dynamics should obey Newton's law. 

(5) The dynamical slide of the block goes on as long as the velocity of the 
block remains positive. When its velocity reaches zero, we assume that the 
block is no more sliding. To account for the heterogeneity and roughness 
of the sliding surface, we assume that the state variable Oq is reset to a 
new random value after the dynamical sliding stops. This random value 
is taken to reflect the characteristics of the new asperities constituting the 
fresh surface of contact. 

(6) After a dynamical slide, the forces exerted by the springs that connect 
the block to its neighbors are updated, as is the new gravitational force 
(if the basal surface has a curvature), the new value of \i is obtained, the 
time counter for frictional creep is reset to zero, and a new process of slow 
frictional creep develops over the new waiting time tf, that is, in general, 
different from the previous one. 

In summary, simulation of the damage process leading to bond rupture between 
blocks proceeds as follows. 

(1) Given an initial configuration of all the blocks within the network, the 
elastic forces exerted by all bonds can be calculated from their exten- 
sion/compression. 

(2) For each bond i subjected to an initial stress so(i), we calculate the cor- 
responding critical time t Cj o(i) at which it would rupture if neither of the 
two blocks connected to it moved in the meantime. For those bonds where 



sq(i) < 5* defined in Equation (3.8), t Cj o(i) is infinite 



(3) Some bonds will eventually fail, modifying the force balance on their blocks 
and accelerating the transition to the sliding regime, after which the stresses 
in the bonds connected to the same blocks are modified. 
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